자 이제 intercept에만 random effects가 주어진 상태에서 어떤 경우에 standard error가 같아지는지 알아보겠습니다. 이를 만족하기 위해서는 다음과 같은 엄격한 조건이 모두 성립해야합니다.
예시 데이터를 확인해보겠습니다. 참고로 예시 데이터의 값들은 모두 임의로 채워넣었기에 결과가 유의하고 않고는 의미가 없음을 미리 말씀드립니다.
이 데이터는 총 20명의 환자를 각 4번씩 반복측정한 데이터입니다. NO열은 각 환자의 ID를 의미하고, Time은 총 4개로 수술 직후, 1개월 후, 6개월 후, 1년 후로 나뉘어 있습니다. Group의 경우 65~79세인 경우(NO가 1-10인 환자)와 80세 이상(NO가 11-20인 환자)인 두 가지 경우가 존재합니다. 이 데이터는 결측이 하나도 존재하지 않고 20명의 환자에 대해 각각 모두 같은 시점에 같은 횟수가 측정된 Complete & Balanced data입니다. 우리는 이러한 반복측정자료에서 Time과 Group에 대한 정보로 Value값에 대해 예측하고 싶어하는 상황입니다. Compounds Symmetry만 만족한다면, 앞서 말했던 것처럼 모든 S.E. 값이 동일하게 나올 것입니다. 한 번 확인해보겠습니다.
lme4패키지의 lmer 함수를 사용하여 진행해보면 아래와 같은 결과를 얻을 수 있습니다.
summary의 Random effects를 보면 intercept에 대해서만 적용된 모습을 확인할 수 있습니다. Fixed effects의 Std.Error을 보게 되면 time effect와 interaction effect에 해당하는 부분의 Std. Error가 모두 같은 모습을 확인할 수 있습니다. 어째서 이런 결과가 나타나는 것일까요?
Linear Mixed Model
지금 우리가 다루고 있는 모델은 여기에 상관관계를 위해 \(u\)라는 term을 추가한 linear mixed model이고, 오차의 독립성 가정을 하지 않기에 선형 회귀 모델과는 조금 다른 구조와 분산-공분산 행렬을 가집니다. 이를 아래와 같이 표현할 수 있습니다.
\(Y = X\beta+Zu+ \epsilon\)
\(u \overset{\mathrm{iid}}{\sim} N(0, \sigma_u^2I), \ \epsilon \overset{\mathrm{iid}}{\sim} N(0, \sigma^2I)\)이고, 이 둘은 독립이라고 가정하겠습니다.
i번째 환자 데이터의 분산-공분산 구조를 살펴보면 다음과 같습니다. \[
Var(Y_i) = V_i = Var(X_i\beta+Z_iu_i+\epsilon_i)
\] \[
= Var(Z_iu_i+\epsilon_i)
\] \[
= Var(Z_iu_i) + Var(\epsilon_i) (\because \epsilon \perp\!\!\!\perp u)
\] \[
= Z_iVar(u_i)Z_i^T + Var(\epsilon_i)
\]
총 20명의 환자가 4번씩 반복측정한 상황이고, random effects가 intercept term에만 존재하기에 \(u_i\)는 univariate normal distribution을 따르는 확률변수가 되고, \[Z_i = \begin{bmatrix}
1 & 1 & 1 & 1
\end{bmatrix} \] 라고 할 수 있으므로 \[
V_i = \begin{bmatrix}
1 & 1 & 1 & 1
\end{bmatrix}
\sigma_u^2
\begin{bmatrix}
1 \\ 1 \\ 1 \\ 1
\end{bmatrix} +
\sigma^2
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix}
\]
\[
= \begin{bmatrix}
\sigma^2+\sigma_u^2 & \sigma_u^2 & \sigma_u^2 & \sigma_u^2 \\
\sigma_u^2 & \sigma^2+\sigma_u^2 & \sigma_u^2 & \sigma_u^2 \\
\sigma_u^2 & \sigma_u^2 & \sigma^2+\sigma_u^2 & \sigma_u^2 \\
\sigma_u^2 & \sigma_u^2 & \sigma_u^2 & \sigma^2+\sigma_u^2
\end{bmatrix}
\] 라는 결과를 얻게 됩니다. 이를 통해 random intercept 모델은 compound symmetry 구조를 만족함을 확인할 수 있습니다.
자 그럼 이제 \(\beta\)의 추정치를 구해보겠습니다. 일반적인 Linear Regression Model에서 OLS(Ordinary Least Square)방법을 이용해 \(\hat\beta=(X^TX)^{-1}X^TY\)를 도출할 수 있다고 했고, 대략적인 과정은 다음과 같습니다.
\[
Y = X\beta + \epsilon
\]
\[
\text{Let }
Q(\beta) = (Y - X\beta)^T(Y-X\beta)
= Y^TY -2Y^TX\beta +\beta X^TX \beta
\] \[
\text{Then } \arg\min_\beta Q(\beta)=\hat\beta \text{ would be the vector that satisfies } \frac{\partial Q(\beta)}{\partial \beta} = 0
\]
\[
\text{since }Q(\beta) \text{ is a convex function w.r.t } \beta.
\]
\[
\text{We know that }X^TX \text{ is symmetric and by assuming it's invertible,}
\]
\[
\frac{\partial Q(\beta)}{\partial \beta} = -2X^TY + 2X^TX\beta = 0
\] \[
X^TX\beta = X^TY
\]
\[
\therefore \hat\beta = (X^TX)^{-1}X^TY
\] 자 이제는 이 방법을 살짝 비틀어서 Linear Mixed Model의 경우를 살펴보겠습니다. 일반적인 형태의 모델은 다음과 같았습니다.
\[
Y = X\beta+Zu+ \epsilon
\]
\[
\text{Let } Zu+\epsilon = \epsilon^*
\] \[
\text{Then }Y = X\beta + \epsilon^*
\] \(Var(Y) = V\)는 symmertic & positive-semi-definite matrix이므로 Cholesky Decomposition에 의해 다음을 만족하는 가역인 행렬 \(\Sigma\)가 존재합니다. \(V = \Sigma \Sigma^T\). 이제 \(\Sigma^{-1}\)을 위 모델의 양변에 곱해봅시다. \[
\Sigma^{-1} Y = \Sigma^{-1}X\beta + \Sigma^{-1}\epsilon^*
\] 이제 아래와 같이 변수들을 새로 정의하면, \[
\Sigma^{-1} Y = \tilde Y,
\] \[
\Sigma^{-1}X\beta = \tilde X\beta
\]
\[
\Sigma^{-1}\epsilon^* = \tilde \epsilon
\]
아래와 같이 OLS의 경우와 아주 유사한 형태의 모델을 새로 정의할 수 있습니다. \[
\tilde Y = \tilde X\beta + \tilde \epsilon
\] 이러한 방법을 OLS를 일반화했다고 하여 Generalized Least Square(GLS)라고 부릅니다. 이제 OLS와 형태가 같기 때문에 \(\hat\beta\)의 추정량은 다음와 같이 도출됩니다.
\[
\hat\beta_{GLS} = (\tilde X^T\tilde X)^{-1}\tilde X^T\tilde Y
\]
\[
= ((\Sigma^{-1}X)^T\Sigma^{-1}X)^{-1}(\Sigma^{-1}X)^T \Sigma^{-1}Y
\]
\[
= (X^T{\Sigma^T}^{-1}\Sigma^{-1}X)^{-1}X^T{\Sigma^T}^{-1}\Sigma^{-1}Y
\]
\[
= (X^T(\Sigma^T\Sigma)^{-1}X)^{-1}X^T(\Sigma^T\Sigma)^{-1}Y
\]
\[
= (X^T(\Sigma\Sigma^T)^{-1}X)^{-1}X^T(\Sigma\Sigma^T)^{-1}Y \ (\because \Sigma^T =\Sigma)
\]
\[
= (X^TV^{-1}X)^{-1}X^TV^{-1}Y
\] 임을 알 수 있습니다. 그렇다면 분산-공분산 행렬은 다음과 같이 유도될 것입니다.
\[
Var(\hat\beta_{GLS}) = Var((X^TV^{-1}X)^{-1}X^TV^{-1}Y)
\]
\[
= (X^TV^{-1}X)^{-1}X^TV^{-1} Var(Y) ((X^TV^{-1}X)^{-1}X^TV^{-1})^T
\]
\[
(X^TV^{-1}X)^{-1}X^TV^{-1}V {V^{-1}}^TX {(X^TV^{-1}X)^{T}}^{-1}
\]
\[
(X^TV^{-1}X)^{-1}(X^TV^{-1}X) (X^TV^{-1}X)^{-1} \ (\because VV^{-1}= I, {V^{-1}}^T=V^T)
\]
\[
= (X^TV^{-1}X)^{-1}
\] 앞에서 보인것 처럼 \(V\)는 Compound Symmertry 구조이고, Complete & Balanced한 데이터를 다루고 있기에 Design matrix인 X의 구조도 변하지 않습니다. 결국 standard error을 결정짓는 것은 방금 전 구한 \(Var(\hat\beta_{GLS})\)의 대각성분과 관련이 있을텐데, 이 값이 같기 때문에 fixed effects의 standard error가 같게 나오는 것입니다.